#include <cstdlib>
#include <cmath>
#include <iostream>
#include <fstream>
using namespace std;
#include "gas_phase_2nd_order.h"
#include "../RateLaw/RateLaw_Kleijn1989.h"

const double TOLERANCE = 1e-15;
const double VERSION = 1.0;

// *******************************************************************
gas_phase_2nd_order::gas_phase_2nd_order(double iKp, double izeta, double itau, double iT, double iP, double it_end, double idt_int, double idt_ext,double iXs_ini, char* iFileName)
{
	int i;
	
	Kp = iKp;
	zeta = izeta;
	tau = itau;
	P = iP;
	T = iT;
	FileName = iFileName;
	
	// Initial value of state variables
	Xs_ini = iXs_ini;
	
	t_begin = 0.0;
	t_end = it_end;
	dt_int = idt_int;
	dt_rec = idt_ext;
	t_node = int(t_end/dt_rec)+1;
	t_arr = new double[t_node];
	for(i = 0;i<t_node;i++)
	{
		t_arr[i] = i*dt_rec;
	}
	
	x[0] = Xs_ini;		// Used in explicit method
	x[1] = 0.0;			

	Xs_sum		 = new double[t_node];
	Xs2_sum		 = new double[t_node];
	DepRate_sum	 = new double[t_node];
	DepRate2_sum = new double[t_node];
	Xi_sum		 = new double[t_node];
	Xi2_sum      = new double[t_node];
	for(i = 0;i<t_node;i++)
	{
		Xs_sum[i]       = 0.0;
		Xs2_sum[i]      = 0.0;
		DepRate_sum[i]  = 0.0;
		DepRate2_sum[i] = 0.0;
		Xi_sum[i]       = 0.0;
		Xi2_sum[i]      = 0.0;
	}

	Repeated = 0;	
	t_idx = 0;
	time = 0.0;	
	Xs = iXs_ini;
	DepRate = ReactionModel(Xs,P,T);

	char FileName2[256];
	strcpy(FileName2,"");
	strncat(FileName2,FileName,strlen(FileName)-4);
	strcat(FileName2,"_detail.dat");
	cout << "gas_phase_2nd_order: opening " << FileName2 << " for output.\n";
	fout_detail.open(FileName2);
	fout_detail << "class: gas_phase_2nd_order\n";
	fout_detail << "Version: gas_phase_2nd_order v" << VERSION << "\n";
	fout_detail << "Created:\n";
	fout_detail << "Kp = " << Kp << "\n";
	fout_detail << "tau = " << tau << "\n";
	fout_detail << "zeta = " << zeta << "\n";
	fout_detail << "Repeated = " << Repeated << "\n";
	fout_detail << "t_begin = " << t_begin << "\n";
	fout_detail << "t_end = " << t_end << "\n";
	fout_detail << "dt_rec = " << dt_rec << "\n";
	fout_detail << "dt_int = " << dt_int << "\n";
	fout_detail << setw(6) << "Run #" << setw(6) << "t(s)" << setw(15) << "Xs";
	fout_detail << setw(15) << "W(layer/s)" << setw(15) << "Xi";
}
// *******************************************************************
gas_phase_2nd_order::~gas_phase_2nd_order()
{
	double Xs_mean;
	double Xs_var;
	double Xs_var2;
	double DepRate_mean;
	double DepRate_var;
	double DepRate_var2;
	double Xi_mean;
	double Xi_var;
	double Xi_var2;
	int i;
	
	fout_detail.close();

	ofstream fout;
	fout.open(FileName);
	fout << "class: gas_phase_2nd_order\n";
	fout << "Version: gas_phase_2nd_order v" << VERSION << "\n";
	fout << "Created:\n";
	fout << "Kp = " << Kp << "\n";
	fout << "tau = " << tau << "\n";
	fout << "zeta = " << zeta << "\n";
	fout << "Repeated = " << Repeated << "\n";
	fout << "t_begin = " << t_begin << "\n";
	fout << "t_end = " << t_end << "\n";
	fout << "dt_rec = " << dt_rec << "\n";
	fout << "dt_int = " << dt_int << "\n";
	fout << setw(6) << "t(s)" << setw(15) << "<Xs>" << setw(15) << "Var(Xs)";
	fout << setw(15) << "<W>(layer/s)" << setw(15) << "Var(W)(layer/s)";
	fout << setw(15) << "<Xi>" << setw(15) << "Var(Xi)" << "\n";
	for(i=0;i<t_node;i++)
	{
		Xs_mean = Xs_sum[i]/Repeated;
		Xs_var2 = (Xs2_sum[i]-Repeated*Xs_mean*Xs_mean)/Repeated;  // Because of numerical error, Xs_var2 CAN be a very negative number sometimes. 
		if(Xs_var2<-1e17)
		{
			cerr << "ERROR@gas_phase_2nd_order::~gas_phase_2nd_order().\n";
			cerr << " i=" << i << "Sum(Xs2)=" << Xs2_sum[i] << "<Xs>=" << Xs_mean << "Var(Xs)^2=" << Xs_var2 << "\n";
		}
		else if(Xs_var2<0 && Xs_var2 >=-1e17)
		{
			Xs_var2 = 0.0;
			Xs_var = 0.0;
		}
		else
			Xs_var = sqrt(Xs_var2);
		
		DepRate_mean = DepRate_sum[i]/Repeated;
		DepRate_var2 = (DepRate2_sum[i]-Repeated*DepRate_mean*DepRate_mean)/Repeated;

		if(DepRate_var2>=0)
		{
			DepRate_var = sqrt(DepRate_var2);
		}
		else if(DepRate_var2<0 && DepRate_var2 >=-TOLERANCE)
		{
			DepRate_var = 0.0;
		}
		else
		{
			cerr << "ERROR@gas_phase_2nd_order::~gas_phase_2nd_order().\n";
			cerr << " i=" << i << "Sum(W2)=" << DepRate2_sum[i] << "<W>=" << DepRate_mean << "Var(W)^2=" << DepRate_var2 << "\n";
			system("pause");
		}

		Xi_mean = Xi_sum[i]/Repeated;
		Xi_var2 = (Xi2_sum[i]-Repeated*Xi_mean*Xi_mean)/Repeated;  // Because of numerical error, Xs_var2 CAN be a very negative number sometimes. 
		if(Xi_var2<-TOLERANCE)
		{
			cerr << "ERROR@gas_phase_2nd_order::~gas_phase_2nd_order().\n";
			cerr << " i=" << i << "Sum(Xi2)=" << Xi2_sum[i] << "<Xi>=" << Xi_mean << "Var(Xi)^2=" << Xi_var2 << "\n";
		}
		else if(Xi_var2<0 && Xi_var2 >=-TOLERANCE)
		{
			Xi_var = 0.0;
		}
		else
			Xi_var = sqrt(Xi_var2);

		fout << setw(6) << setprecision(6) << t_arr[i];
		fout << setw(15) << setprecision(6) << Xs_mean;
		fout << setw(15) << setprecision(6) << Xs_var;
		fout << setw(15) << setprecision(6) << DepRate_mean;
		fout << setw(15) << setprecision(6) << DepRate_var;
		fout << setw(15) << setprecision(6) << Xi_mean;
		fout << setw(15) << setprecision(6) << Xi_var;
		if(i<t_node-1 )
		{
			fout << "\n";
		}
	}
	fout.close();

	delete [] t_arr;
	delete [] Xs_sum;
	delete [] Xs2_sum;
	delete [] DepRate_sum;
	delete [] DepRate2_sum;
	delete [] Xi_sum;
	delete [] Xi2_sum;
}
// *******************************************************************
double gas_phase_2nd_order::update_explicit(double Xin)
{
	if(Xin > 1.0||Xin < 0.0) {cerr << "ERROR@gas_phase_2nd_order::update_explicit: Invalid input. Xin = " << Xin << "\n";abort();}
	
	// Use explicit Euler
	x[0] = x[0]+dt_int*x[1];
	x[1] = x[1]+dt_int*(Kp*Xin-2*zeta*tau*x[1]-x[0])/(tau*tau);
	Xs = x[0];
	if (Xs < 0.0 || Xs > 1.0) {cerr << "ERROR@gas_phase_2nd_order::update_explicit: Invalid Xs. Xs = " << Xs << "\n";abort();}
	
	DepRate = ReactionModel(Xs,P,T);
	time += dt_int;
	
	if (t_idx >= t_node) {cerr << "t_idx exceed bounary.";abort();}
	if(abs((time-t_arr[t_idx])/dt_rec)<0.1)
	{
		Xs_sum[t_idx] += Xs;
		Xs2_sum[t_idx] += Xs*Xs;
		DepRate_sum[t_idx] += DepRate;
		DepRate2_sum[t_idx] += DepRate*DepRate;
		Xi_sum[t_idx] += Xin;
		Xi2_sum[t_idx] += Xin*Xin;
		
		fout_detail << "\n";		
		fout_detail << setw(6)  << setprecision(6) << Repeated;
		fout_detail << setw(6)  << setprecision(6) << time;
		fout_detail << setw(15) << setprecision(6) << Xs;
		fout_detail << setw(15) << setprecision(6) << DepRate;
		fout_detail << setw(15) << setprecision(6) << Xin;

		t_idx++;
	}
	return DepRate;
}
// *******************************************************************
void gas_phase_2nd_order::get()
{
	cout << "class: gas_phase_2nd_order\n";
	cout << "Created:\n";
	cout << "Kp = " << Kp << "\n";
	cout << "tau = " << tau << "\n";
	cout << "zeta = " << zeta << "\n";
	cout << "Repeated = " << Repeated << "\n";
	cout << "t_begin = " << t_begin << "\n";
	cout << "t_end = " << t_end << "\n";
	cout << "dt_rec = " << dt_rec << "\n";
	cout << "dt_int = " << dt_int << "\n";
}
// *******************************************************************
void gas_phase_2nd_order::reset()
{
	Repeated++;
	time = 0.0;
	t_idx = 0;
	Xs = Xs_ini;
	DepRate = ReactionModel(Xs,P,T);

	// Gangshi: x[0] and x[1] not been reset here, this might be the reason, since Xs is later updated by x[1]
	x[0] = Xs_ini;		// Used in explicit method
	x[1] = 0.0;	

	Xs_sum[t_idx]		+= Xs;
	Xs2_sum[t_idx]		+= Xs*Xs;
	DepRate_sum[t_idx]	+= DepRate;
	DepRate2_sum[t_idx] += DepRate*DepRate;
	Xi_sum[t_idx]		+= 0.0;
	Xi2_sum[t_idx]		+= 0.0;
	t_idx++;
}
// *******************************************************************
int gas_phase_2nd_order::notfinish()
{
	return t_idx < t_node;
}
